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Qh! Abstract 

D , We numerically study the effects of high gluon density at small x on the evolution 



of gluon distribution function in both hadrons and nuclei. Using a newly derived, 
Wilson renormalization group-based evolution equation which includes n to 1 gluon 



5h ' ladder fusion, we find significant reduction in nuclear gluon distribution function 

for large nuclei at zero impact parameter at energies relevant for RHIC and LHC 
experiments. 



1 Introduction 

Gluons are the most abundant partons in hadrons and nuclei at small x (high energy) and 
as such, will determine the behavior of many physical observables such as hadronic/nuclear 
cross sections through their initial distribution. Once the initial distribution of gluons in a 
hadron or nucleus is known, its change with x and Q^ can be predicted using the powerful 
machinery of perturbative QCD through the standard QCD evolution equations such as 
DGLAP ^ and BFKL §]. Both DGLAP and BFKL evolution equations can describe the 
available experimental data quite well in a fairly broad range of x and Q^ with appropriate 
parameterizations. It is interesting to note that both equations predict a sharp growth 
of the gluon distribution function as x grows smaller which is clearly seen in the Deep 
Inelastic Scattering (DIS) experiments at HERA (see [^ for a recent review). 

This sharp growth of the gluon distribution function will have to eventually slow down 
in order to not violate unitarity (Froissart [Q) bound on physical cross sections. Gluon 
recombination is believed to provide the mechanism which is responsible for this slow down 
or a possible saturation of the gluon distribution function at small x. In other words, the 
number of gluons at small x will be so large that they will spatially overlap and therefore, 
gluon recombination will be as important as gluon splitting and the standard evolution 
equations like DGLAP will have to be modified in order to take this into effect. A first 
step in this direction was taken in ^ by Gribov, Levin and Ryskin (GLR) who suggested 
the form of the first non-linear correction to DGLAP and identified the diagrams which 
contribute. In 0, Mueller and Qiu (MQ) made a Gaussian like ansatz for the gluon 4- 
point function which was taken to be proportional to the square of the gluon distribution 
function (2-point function). They then proceeded to calculate the numerical coefficient of 
the non-linear term. 

It was shown in ^ that assuming reasonable values for R, a phenomenological pa- 
rameter which could be taken to be either the proton or valence quark radius, gluon 
recombination effects were negligible in hadrons for not too small values of x but could 
be significant for large nuclei. It is perhaps helpful to mention that GLR/MQ and gluon 



recombination based approaches in general are formulated in the infinite momentum frame. 
There has been much work inspired by the approach of GLR/MQ which show that gluon 
recombination leads to saturation of gluon density at small x 0]. However, GLR/MQ 
approach includes only the first non-linear term in the evolution equation and will not be 
valid at very small x where contribution of higher order terms will be as important as the 
first order correction and one will need to include them as well. Also, in the original ap- 
proach of GLR/MQ there is no information about the impact parameter dependent gluon 
distributions and one typically assumes a factorization of the impact parameter and the 
usual gluon distribution function at all x. 

In p, 1^, a new evolution equation for gluon distribution function (or any gluonic n- 
point function) was derived which is valid in the small x region. It was shown in |]lOl that 
the new equation reduces to all the standard evolution equations like BFKL, DLA DGLAP 



and GLR/MQ in the low density limit. In ITTl double leading log limit of the new equation 
was considered and a closed form was obtained which generalizes the GLR/MQ equation. 
It was shown that this new equation slows down the growth of gluon distribution function 
at small x consistent with unitarity limits. In this paper, we will numerically solve this 
equation and investigate high gluon density effects on the evolution of gluon distribution. 
There are many interesting situations where understanding these effects should be useful. 

Mini-jet production at high energy is an example where high gluon densities will play 
an important role. Mini-jets will be important at RHIC and will dominate at LHC over 
soft phenomena. Nuclear shadowing of initial gluon distribution and high gluon density 
could significantly reduce the initial mini-jet and total transverse energy production. Such 
reduced initial energy density will also affect the subsequent parton thermalization. 

Another example is heavy quark production where high gluon density effects may 
make a dramatic difference specially at LHC. Since the probability for making a heavy 
quark pair is proportional to square of gluon density, any depletion in number of gluons 
will make a significant difference in the number of heavy quark pairs produced. 

In section 2, we discuss the relation between gluon distributions in nuclei and hadrons 
followed by a brief description of Wilson renormalization group and effective action ap- 
proach to high density/small x QCD. In section 4, We outline the semi-classical approach 



for solving the general equation and use numerical methods to solve them. We finish by a 
discussion of our results and their experimental implications as well as the limitations of 
our approach. 

2 Gluons in Hadrons and Nuclei 



Gluon distribution function in a hadron, ^(^(x, Q^), has been studied quite extensively. 
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Theoretically, once the initial distribution function at a given scale Xq and Qq is known 



one can calculate the distribution function at a different scale x and Q^ using the standard 
perturbative QCD-based evolution equations, for instance, the DGLAP equation. How- 
ever, the initial distribution is non-perturbative and has to be supplied as an input to the 
evolution equation and is usually taken from parameterized experimental data. Alterna- 
tively, one can use the BFKL equation to study evolution of gluon distribution function 
with X. It is interesting to notice that both of these evolution equations predict a sharp 
rise of the gluon distribution function which would eventually lead to violation of unitarity 
(Froissart bound @]). 

Gluon distribution function of a proton can be measured indirectly in DIS experiments 
at HERA and elsewhere by measuring virtual photon-proton cross section a^*p where 

^rp^ -^^^i^^Q^)- (1) 

It should be kept in mind that eq. ([^) is a leading twist relation and will break down when 
we consider higher twist terms in the evolution of the gluon distribution function. The 
theoretically predicted sharp growth of the gluon distribution function with x is observed 
experimentally at all Q^ as well as at fixed (and small ) x with increasing Q^. Even though 
DGLAP evolution equation fits the data very well, one can also use BFKL to explain the 
data and as of now, one can not experimentally distinguish between the two scenarios. 

This sharp growth of gluon distribution function is expected to slow down eventually 
due to mutual interactions between gluons when they start to spatially overlap. This is 
usually referred to as saturation of gluon density and will happen when probability of two 



gluons recombining into one is as large as the probability for a gluon to split into two 
gluons. In other words, when 

«., xG ^ ^2) 
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one has to include recombination effects which are neglected in DGLAP and BFKL evolu- 
tion equations. 

To illustrate this, one may write a generic evolution equation for gluon distribution 
function in the high density region as 
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where y = Inl/x and ^ = InQ^. The first term of the sum on the right hand side of 
eq.(^) is just the DGLAP equation whereas the second term is referred to as GLR/MQ 
since it was first investigated in [^ and its numerical coefficient was calculated by Mueller 
and Qiu |Q. An exact and formal evolution equation for gluon distribution function to 
all orders in gluon density in the high density (small x) region was first derived in [§ in 
the infinite momentum frame. Whereas BFKL equation can be thought of as a ladder of 
reggized gluons, the general equation derived in can be thought of as having n ladders of 
reggized gluon fusing into one and in this sense, it is the generalization of BFKL equation 
appropriate for the high gluon density region. 

Until recently, there was no evidence that high density effects in a hadron were experi- 
mentally relevant even at the smallest x (highest energies) achieved at HERA in agreement 
with estimates made in MQ if one assumed reasonable values for the hadron/quark radius 
R. It should be emphasized that R is a purely phenomenological parameter and can not 
be derived from perturbative QCD. There is a very recent report on slope of the proton 
structure function Ff from HERA ||12| which may be an indication of importance of higher 



twist effects in proton |T3 



Gluon distribution function in a nucleus is intimately related to that in a hadron. 
Typically, one assumes that nucleus is a weakly bound system of nucleons so that one can 
neglect inter nucleon forces which is equivalent to taking the nucleus to be a dilute system 



of nucleons in the transverse plane. This is a good approximation for hard processes in 
high energy nuclear collisions under normal conditions. Furthermore, if one assumes that 
density of gluons in a hadron is low, then one can simply relate distribution of gluons in 
nuclei and hadrons by 

xG^{x,Q^) = AxG{x,Q^) (4) 

where A is the atomic mass number, xG^{x, Q"^) and xG{x, Q^) are the nucleus and hadron 
(proton) gluon distribution functions. 

Even if high density effects are not well established in hadrons, they are expected to be 
much more important for heavy nuclei. Non-linear terms in the evolution equation for gluon 
distribution function in a nucleus become appreciable at a larger x (lower energy) than for 
hadrons. In this sense, nucleus can be thought of as an amplifier of non-linear effects in 
QCD. These high density effects may very well be present in experiments planned at RHIC 
and LHC which underscores the crucial importance of a theoretically well-defined approach 
to nuclear gluon distribution function. Also, having nuclear beams at HERA would be of 
great help pinning down these effects and would be complementary to experiments planned 
for RHIC and LHC. 

One of the advantages of our approach is that it can be used to investigate the gluon 
distribution function in both hadrons and nuclei and its impact parameter dependence 
without any assumptions on the form of impact parameter at small x. This will allow a 
systematic and rigorous determination of the change in the impact parameter of nuclear 
gluon distribution function with energy. To our knowledge, this is the first time that x 
and Q^ dependence of nuclear gluon distribution function as well as its impact parameter 
dependence have been derived from QCD in the high density region. 

3 The General Evolution Equation 

In this section we will briefly review the Wilson renormalization group and effective action 
approach to small x QCD as developed in [Q-|]TT|, |T^-[0. To make this paper self- 



contained, we will include some of the results reported above as needed. A few years ago, 
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McLerran and Venugopalan WM suggested that for very large nuclei and/or at very small x 



one can use weak coupling, semi-classical methods to calculate structure functions. They 
considered a large nucleus in the infinite momentum frame and argued that as long as 
number of valence quarks per unit area per unit rapidity is large, they can be treated as 
static, classical sources of color charge to which the long wavelength gluonic fluctuations 
(small X gluons) couple. In order to perform color averaging over the hadron/nucleus state, 
they assumed a Gaussian weight for color configurations. They proceeded to solve Yang- 
Mills equations of motion and calculated the gluon distribution function in lowest order 
in the coupling constant. Quantum corrections to the classical result were computed in 
T5| and analogous to standard perturbation theory, large logarithmic factors (In l/x) were 



encountered which necessitated a formalism which would resum these large logarithmic 
factors in presence of a non-trivial background (classical) field. 

In [0], McLerran- Venugopalan action was generalized as the following 



S = -^^j d'xG>:^G% + i j d\^F[p\x^)] 

+ — (fx_Ldx' 6{x')p''{x_i_)tTTaW_oc,oc,[A']{x' , X±) 

where W is the Wilson line in the adjoint representation along the x^ axis 



W-oo,oo[-A ]{x ,Xt) = Pexp 



ig / dx^A^{x^,x ,Xt)Ta 



(5) 



The nucleus/hadron is represented by an ensemble of color charges localized in the plane 

x^ = with the (integrated across x^) color charge density p{xj_). Statistical weight of a 
configuration p{x±) is given by 

Z = exp{-F[p]} (6) 

In light cone gauge A'^ = and at the tree level, the chromoelectric field is determined by 
the color charge density through the equations 

G+' = -6{x-)ai{x^) (7) 
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where the two dimensional vector potential ai{x±) is "pure gauge" and is related to the 
color charge density by 

9i«; - dja^ - r'^'a'la'j = 

d^at = -P^ (8) 

One can then consider quantum fluctuations in background of this classical field and 
separate hard and soft modes (in light cone longitudinal momenta) of the fluctuations, 
keeping terms quadratic in hard fluctuations. Integrating out the hard modes generates 
the renormalization group equation which has the form of the evolution equation for the 

statistical weight Z p 

In the notation used in Eq. @, both u and v stand for pairs of color index and trans- 
verse coordinate with summation and integration over repeated occurrences implied. The 
evolution in this equation is with respect to the rapidity y, related to the Bjorken x by 
y = In 1/x. 

The quantities x[p] and a[p] have the meaning of the mean fluctuation and the average 
value of the extra charge density induced by the hard modes of quantum fluctuations. They 
are functionals of the external charge density p. The explicit expressions have been given 
ini. 

Using this equation for the statistical weight Z, one can derive evolution equations for 



n-point functions of gluon field M. In [11|, double leading log limit of the 2-point function 



(gluon distribution function) was investigated and shown to be 

- aa 

(10) 



4- < atiX)a^{Y) >= 4a, 
dy 



< ^'larai^' > 



In the high density limit where a^ ^ 9^, one can neglect the derivative term in 
the denominator above and the right hand side is a constant which leads to the gluon 
distribution function growing only logarithmically with x (energy) consistent with unitarity. 
In the low density limit where a^ ^ d'j_, one can expand the denominator in the above 

8 



equation. The first term of tlie expansion gives tlie DGLAP equation. Furthermore, if one 
assumes a factorization of the 4-point function in terms of the 2-point function (as assumed 
by GLR/MQ), one recovers the GLR/MQ equation [Q. This is equivalent to ignoring all 
correlations between gluon fields except that they are constrained to be in an area of ttR^. 
With these assumptions, one can actually perform the color averaging in ([T0|) which leads 
to 
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dyd^ 



xG{x, Q, 6j 



NcjN, - 1) ^2 



l--exp(-)Ei(-; 

Kj Kj Kj 



:ii) 



where 



2a. 



;xg{x,Q,b±) 



7r{N, - l)g2 
and Ei(x) is the exponential integral function defined as ]T7[ 

-(l+i)x 



(12) 



^00 g (.J-T^i-;-^ 

Ei(x) = dt , X > 

Jo 1 + 1 



(13) 



In the low density limit, one can expand equation (|TT|). Keeping the first term, we get 



xG{x, Q, b±) = — ^—^ xG{x, Q, bj 



dydi 



TC 



(14) 



which is the DLA DGLAP (small x limit of DGLAP) equation. In the high density limit. 



eq. (|TT|) gives 



d^ 



xG{x, Q, b_[ 



N,{N,-l)^^ 



which leads to a gluon distribution of the form 



-Q' 



(15) 



xG{x, Q, b±) ~ Q^ In 1/x. 



(16) 



Let's consider the impact parameter dependent gluon distribution function 
xG{x, Q, b±) which is related to gluon distribution function xG{x, Q) by 



xG{x,Q) = / (fb± xG{x,Q,b±). 



(17) 



It is usual to factor out the impact parameter dependence of the distribution function and 
write xG{x, Q, b±) = S{b±) xG{x, Q) where 5'(6j_) is the nucleus/nucleon shape function 
and can be taken to be a Gaussian 



,2 /r2 



SC-x) = ^j^ (18) 

SO that J d'^b± S{b±) = 1. This factorization introduces the phenomenological parameter 
R which is taken to be the nuclear /hadronic radius. As long as one is using the DGLAP 
evolution equation, this parameter does not come into play since DGLAP equation is 
linear in gluon density. However, as we consider the first non-linear term in the evolution 
equation as in GLR/MQ equation, parameter R becomes relevant and one needs to define 
it precisely. Since in GLR/MQ impact parameter dependence is factorized, one can only 
make plausible estimates of R. 

In general, this factorization of impact parameter will break down with evolution in 
X and Q^ simply because gluon densities are expected to be higher in the central {bj_ = 0) 
region than the peripheral {b± ~ R) region and so therefore will evolve differently. This 
would lead, in the general case, to a breakdown of factorization of impact parameter and 
Gaussian ansatz for the nucleus/nucleon shape function. In the present case where we 
are working in the double logarithmic region, the Gaussian ansatz for the shape function 
should still hold but the area (or radius R) would change with x and Q^. This basically 
amounts to the rise of perturbative cross sections with energy. In this work, we factorize 
the impact parameter only at the starting point of our evolution xq and Qq where non- 
linear effects are believed to be experimentally absent. The evolution equation will then 
predict the change of this "area" with energy. 

4 Solving the General Equation 

In this section we will outline the procedure to numerically solve the general equation ([Til). 
We will use the method of characteristics which converts a partial differential equation to 
a set of coupled ordinary differential equations [l^ (see also [^ |19|, |2^ for an illustration of 
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this method). We will also use MQ [Q normalization of 4-point function in terms of 2-point 
functions in order to facilitate comparison of our results with those where one includes only 
the first non-linear term in the evolution equation. This amounts to a simple rescaling of 
our gluon distribution function 



xG{x, Q, 6j 



6 



xG{x, Q, 6j 



(19) 



In the following, we will closely follow the derivation of Ayala et al. in |0 and rewrite (|TT 
in terms of the density factor k so that one gets 

Ncas 



dyd^ dy 

where the rescaled k is now 



vr 



Nra, n^ 



K 



1 - -Exp{-)Ei{-) 

K K K 



xG{x, Q, 6j 



(20) 



TT 3g2- 

In the semi-classical approximation, one can write the solution to (pO|) as 



(21) 



K = e 



(22) 



and neglect 



d^S dS dS 



dyd$, dy d^ 



Defining ^ = u and H =7, we get 



9€ 



u;(7 + l) = $(5) 



(23) 



where 



HS) 



7r 



1-e 



''^'-' E,(e-') 



In terms of these variables, the set of characteristic equations become 

1 , d'j 79$ 



dS 27 + 1 d^ 



dy (7 + 1)^ 



dy (7 + 1)2 






dy -f + ldS 



(24) 



(25) 
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Notice that these equations are identical in form to those in |^ except that our function 
$(5*) is different and will therefore result in different solutions. 

In order to solve these equations, We need some initial conditions. Since they are first 
order ordinary differential equations, we will need to specify their initial values denoted 
by 5*0, 7o and ,^0 at some initial point yo- In order to clarify these initial conditions, it is 
helpful to write them explicitly in terms of the gluon distribution function 



So = In 



vr 2Q 



xG{xo,Qo,bj 



d 
70 = ^ ^^^(a;,^,^^)!^^,^^ -1 

eo = Ing^ (26) 

We will choose the initial Xq and Qq such that the non-linear terms are negligible in the 
evolution equation. For a nucleon, this requirement is not very restrictive since non- 
linear effects are small in a broad range of x and Q^. In a nucleus, however, it is known 
experimentally that there is a narrow range of x such that the shadowing ratio S = -^ ~ 1 
so that we will restrict our initial point Xq and to lie in this region. From experimental 



data |21|, ^, it appears that xq ~ 0.05 — 0.07 is a reasonable value so that for the sake 
of definiteness, we will choose xq = 0.06 but have checked that our results are not very 
sensitive to variation of Xq in this range. We also choose Qo = 0.7 (in practice, most 
characteristic lines start at a higher Qq) for the following reasons: quite surprisingly, all 
HERA data can be explained by starting at such low values of Qq so that perturbative 
QCD seems to hold at such small values in DIS (GRV parameterization of parton densities 
start at Qq = O.bGeV). Also, since we are using the method of characteristics to solve these 
equations, it is useful to start from a low Qq in order to be able to find the characteristic 
lines in a broad range of x. Most importantly, we want to get an upper limit on the amount 
of perturbative shadowing generated so that it is helpful to start from as low virtualities 
as allowed by perturbative QCD. 

Having chosen our initial point Xq and Qq, we use CTEQ parameterization of the 
proton gluon density to get XqG^^Xq, Qq). Also, at the initial point Xq and Qq, we assume 
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that factorization of the impact parameter is vahd as discussed in some length earher 

xoG{xo, Qo, b±) = S{b±)xoG{xo, Qo). 

Putting everything together, at the initial point xq and Qo, the impact parameter 
dependent gluon distribution function in nuclei and hadrons can be written as 

xoG^ixo, Qo, h^) = A^-^-^xoC'ixo, Qo) (27) 

and 

xoG^lxo, go, &±) = 5^xoG^(xo, go) (28) 

TTiX 

where R\ and R^ are the nuclear and hadronic areas at the initial point xo and go- This 
completely fixes our initial conditions for solving the set of coupled ordinary differential 
equations in (^). 

We would like to emphasis that the factorization of the impact parameter dependent 
distribution into a Gaussian shape function S{h^_) and the standard gluon distribution 
function xG{x, Q) is done only at the initial point and in principle will not hold when one 
goes to small x where solution of the evolution equation will determine its functional form. 
Here, we will mostly work at zero impact parameter since that is where gluon densities and 
hence non-linearities are most important but we intend to investigate impact parameter 
dependence of our results in more detail in future work. 

Having determined our initial conditions, we use the 4th order Runge-Kutta method to 
solve the set of characteristic equations in (P^). Some of the characteristic lines are shown 



in Figure (1) for illustration. All the lines start at Xq = 0.06 and end at g = lOGeV. In 
order to find the gluon distribution function xG{x, Q,b±) at a given x, one would need to 
vary the initial go until the corresponding characteristic passing through the given x and Q 
is found. For the range of x and Q considered here, variation of go is between 0.8 — l.bGeV 
so that even though the evolution formally starts at a low value of go = 0.7GeV, the actual 
initial go is higher. 
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Figure 1: Some characteristic lines of equation ( pSf ) starting at Xq = 0.06 (yo = 2.81 



In Figures (2) and (3), we show the ratio 



R{x,Q,b_i 



(29) 



xGDG^^P{x,Q,b^) 

at b± = for both A = I and A = 200 at Q = 2 GeV and Q = 5 GeV. We have also 
taken Ra = 5 fm, R = 1 fm and as = .25. Here, xG"'^^^ refers to solution of equation 
(p!l|) while xG^^^"^^ is the solution of (Q). For a proton, we get a 15 — 20% reduction in 
the number of gluons at x ~ 10~^ as compared to DLA DGLAP while for a Gold or Lead 
nucleus, there is a 50 — 55% reduction at a; ~ lO^''. It is expected that these results will 
have some dependence on the numerical values of the hadron or nucleus radius as well as 
the coupling constant a^. For example, one can find values like R = .bfm for radius of 
proton and so on in the literature. Also, one could use a more realistic shape function like 
Woods-Saxon rather than a Gaussian but these changes are not expected to make more 
than a few percent change in our results. In this work, we are interested in the overall 
features of the non-linear effects and will investigate these details later. 

In Figure 4, we show the Q dependence of this ratio at different x for A = 200. To 
make comparison easier, we normalize R = 1 a.t Q = 1. For a proton, the Q dependence 
was found to be negligible and is not shown. 
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Figure 2: R{x, Q, b±), as defined in (|29|) at Q = 2 GeV and b± = 0. 




o 0.6 i- 

II 

in- 0.5 I- 

II 

a 

£- 0.4 I 



0.3 r 
0.2 i- 
0.1 i- 



0.0 
0.0001 



A=1 
A=200 



0.0010 



0.0100 



0.1000 



Figure 3: Same as in Figure (|) at Q = 5 GeV. 
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Figure 4: R[x, Q, b±) as a function of Q at different x for A = 200. 



As is seen, the non-linearities become less important at higher Q. However, the low 
Q behavior is peculiar since one expects a monotonous decrease of R with decreasing Q 
while it is seen to eventually increase with decreasing Q and tend to 1. This dip in Q is 
a consequence of using the method of characteristics to solve partial differential equations 
and is not an artifact of our formalism (see, for example. Figs. 23, 24 in 0). 

To find the gluon distribution function xG{x, Q) at Q, one needs to find the charac- 
teristic line of the corresponding evolution equation passing through that value of Q by 
finding the value of Qq from which the desired characteristic line starts. xG"^^^^ and 
^qdglap ga^^igfy different evolution equations and therefore have different characteristic 
lines which start at different values of Qo in order to reach the point Q. In other words, 
the initial starting Qq is never exactly the same for the two evolution equations. This 
would not be important if the distribution functions at the initial Qo were slowly varying 
which is not the case for xG'^^^^{x,Q) since it includes high density effects. To make a 
completely self-consistent treatment of perturbative shadowing possible, one would need to 
start from parameterizations of parton distributions which include initial non-perturbative 
shadowing such as [|^ rather than CETEQ, GRV or MRS. We intend to do this when we 



16 



study perturbative gluon shadowing in nuclei in the near future. 

In Figures (5) and (6), we show the A and b± dependence of our results at fixed Q for 
different values of x. As is seen, non-linear effects set in rather quickly, specially at low x, 
after which they increase slowly. 




20 40 60 80 100 120 140 160 180 200 
A 



Figure 5: R{x, Q, b±) as a function oi A at Q = 5GeV for different x. 



These figures clearly show the importance of the non-linear terms, specially for a large 
nucleus, at values of x which will be reached in the upcoming experiments at RHIC and 
LHC. These non-linear effects will be manifest in terms of shadowing of nuclear gluon 
distribution function and will have to be taken into account at future high energy heavy 
ion experiments. We are currently investigating shadowing of nuclear gluon distributions 

using this formalism. 



5 Discussion 

We have investigated the x and Q dependence of gluon recombination effects on the evolu- 
tion of impact parameter dependent gluon distribution function in hadrons and nuclei using 
a general evolution equation which takes n —>■ 1 gluon ladder fusion into account. Using the 
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Figure 6: R{x,Q,b±) as a function of b± at different values of x. 



method of characteristics, we numerically solved the general evolution equation and found 
that gluon recombination effects are very important specially for central ultra-relativistic 
heavy ion collisions coming up at RHIC and LHC 

A detailed knowledge of nuclear gluon distribution will be essential in understanding 
the outcome of RHIC and LHC experiments. With our formalism and with the solution 
to the general evolution equation at hand, we can investigate several aspects of these 
distributions. The first thing to consider is shadowing of nuclear gluon distributions using 



the present formalism. In addition to the usual shadowing ratio defined as 



xG^{x,Q) 
AxG{x,Q) ' 



we can 



investigate impact parameter dependence of this ratio. This will give information about 
shadowing of gluon distributions in the peripheral as well as the central region. Studying x 
and Q dependence of the average impact parameter will determine how the effective area 
(cross section) of a nucleus changes with increasing energy of nuclear collisions. The A 
dependence of shadowing ratio at small x can also be determined rigorously without any 
model dependent assumptions 0. One should also integrate over the impact parameter 
eventually in order to make possible a comparison with other approaches which take only 



the first non-linear term into account. Work in this direction is in progress ||24 
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One could also consider the effect of shadowing of nuclear gluon distributions on 



physical cross sections like Drell-Yan |^ and heavy quark production. To do this, one 
would need to include the relevant higher twist effects not only in the evolution of gluon 
distribution (as is done here) but also in the relation between the cross section (or F2) and 
gluon distribution function. In |^, the all twist structure function, F2, was computed in 
the infinite momentum frame (see also ||^, ^ for a similar calculation in the lab frame) so 
that all one has to do is to merge the two results [ \i8\ \ . 

There are a few issues which need to be analyzed further. First, one should investigate 
our choice of the initial virtuality Qo in solving the evolution equation. We took a low value 
of Qo = 0.7 GeV (see the comments before eq. (|29D) because it seems to be consistent 
with DIS data from HERA and in order to get maximum perturbative shadowing of gluon 
distribution possible. However, from Figure (4), we see that R decreases rather sharply at 
low values of Q before it starts to increase with increasing Q as it must as a high twist 
effect. This means that our treatment, strictly speaking, is not self-consistent at low values 
of Q. In other words, starting from high values of Q and decreasing Q, we should have a 
monotonous decrease of R as is seen in Figure (4). Further decrease of Q leads to R in- 
creasing in order to match our initial condition that R{xo, Qo) = 1- This indicates that one 
should include some initial non-perturbative shadowing so that R{xo, Qo) 7^ 1. Unfortu- 
nately, this requires detailed knowledge of x, Q, A and b± dependence of gluon distribution 
function at the initial point which necessitates use of many model dependent assumptions. 
There are a number of approaches including vector meson dominance, Pomeron exchange 
models, etc. ( see p9| and references therein) which one could adopt to address this issue. 
This is beyond the scope of present work and will be pursued in future. 

We are also going to study the dependence of our results on our choice of CTEQ 
parameterization of gluon densities by repeating this calculation using other available pa- 
rameterizations like MRS and GRV. However, since we used the parameterized gluon dis- 
tributions at a fairly high initial Xq ~ 0.06 and because parameterization dependence of 
gluon distribution function becomes noticeable only at small values of x, we do not believe 
our results will be sensitive to choice of parameterization. 
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Another point to keep in mind is that we have been working in the leading log ap- 
proximation and therefore have taken a^ to be a constant. It is known that DLA DGLAP 
with fixed as overestimates the gluon number density at small x and to get good agreement 
with experimental data from HERA, one needs to include next to leading order corrections 
to DGLAP. However, our general evolution equation is derived in the leading log approxi- 
mation and will be modified if one goes beyond the leading log approximation. One such 
modification due to next to leading order corrections is to cause running of a^, but this may 
not be the only effect or even the dominant effect so that to be theoretically consistent, 
we have kept everything at the leading log approximation level. However, so long as we 
are working with ratios of distributions, we think next to leading order corrections to this 
ratio will not be large and a leading order calculation such as this should be adequate. 

Our emphasis in this work has been on theoretical self-consistency so that we have not 
made any attempt to fit or reproduce experimental data. The main reason is that until now, 
there has been little effort made to calculate nuclear gluon distribution function directly 
from QCD without resorting to elaborate modeling. Also, there is not much experimental 
data available on nuclear gluon distribution function. What is experimentally measured is 
the structure function F2 and to get the gluon distribution function, one takes logarithmic 
derivative of F2 



d 



^G{x,Q)r~.^^——F2. (30) 



As mentioned earlier, this is a leading twist relation and will not be valid in a general 
all twist calculation such as ours. By using eq. (|30|) to extract the gluon distribution 
function experimentally, one is implicitly assuming that higher twist effects are not present. 
To be theoretically consistent, one should figure out how to extract gluon densities from 
experimental data allowing for higher twist effects. Until this is done, in our opinion, 
one should develop a self-consistent approach derived from fundamental theory with well 
defined approximations and with as less model dependence as possible. This is the approach 
adopted in this work. 
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